Importation de la base:

library(readxl)
Book1 <- read_excel("C:/Users/sanat/OneDrive/Bureau/Hello/Atelier/Book1.xlsx")
BD <- Book1

Chargement des librairies:

library(fpp3)
## ── Attaching packages ────────────────────────────────────────────── fpp3 0.5 ──
## ✔ tibble      3.2.0     ✔ tsibble     1.1.3
## ✔ dplyr       1.1.0     ✔ tsibbledata 0.4.1
## ✔ tidyr       1.3.0     ✔ feasts      0.3.1
## ✔ lubridate   1.9.2     ✔ fable       0.3.3
## ✔ ggplot2     3.4.2     ✔ fabletools  0.3.2
## ── Conflicts ───────────────────────────────────────────────── fpp3_conflicts ──
## ✖ lubridate::date()    masks base::date()
## ✖ dplyr::filter()      masks stats::filter()
## ✖ tsibble::intersect() masks base::intersect()
## ✖ tsibble::interval()  masks lubridate::interval()
## ✖ dplyr::lag()         masks stats::lag()
## ✖ tsibble::setdiff()   masks base::setdiff()
## ✖ tsibble::union()     masks base::union()
library(zoo)
## 
## Attachement du package : 'zoo'
## L'objet suivant est masqué depuis 'package:tsibble':
## 
##     index
## Les objets suivants sont masqués depuis 'package:base':
## 
##     as.Date, as.Date.numeric
library(forecast)
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
## 
## Attachement du package : 'forecast'
## Les objets suivants sont masqués depuis 'package:fabletools':
## 
##     accuracy, forecast
library(tseries)
library(astsa)
## 
## Attachement du package : 'astsa'
## L'objet suivant est masqué depuis 'package:forecast':
## 
##     gas
library(dplyr)
library(lubridate)
library(tsibble)
library(tsibble)
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats 1.0.0     ✔ readr   2.1.4
## ✔ purrr   1.0.1     ✔ stringr 1.5.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter()     masks stats::filter()
## ✖ tsibble::interval() masks lubridate::interval()
## ✖ dplyr::lag()        masks stats::lag()
## ℹ Use the ]8;;http://conflicted.r-lib.org/conflicted package]8;; to force all conflicts to become errors
library(ggplot2)
library(gridExtra)
## 
## Attachement du package : 'gridExtra'
## 
## L'objet suivant est masqué depuis 'package:dplyr':
## 
##     combine
library(grid)
library(ggthemes)
library(corrplot)
## corrplot 0.92 loaded
library(lmtest)
library(sandwich)
library(TSA)
## Registered S3 methods overwritten by 'TSA':
##   method       from    
##   fitted.Arima forecast
##   plot.Arima   forecast
## 
## Attachement du package : 'TSA'
## 
## L'objet suivant est masqué depuis 'package:readr':
## 
##     spec
## 
## Les objets suivants sont masqués depuis 'package:stats':
## 
##     acf, arima
## 
## L'objet suivant est masqué depuis 'package:utils':
## 
##     tar
library(car)
## Le chargement a nécessité le package : carData
## 
## Attachement du package : 'car'
## 
## L'objet suivant est masqué depuis 'package:purrr':
## 
##     some
## 
## L'objet suivant est masqué depuis 'package:dplyr':
## 
##     recode
library(aTSA)
## 
## Attachement du package : 'aTSA'
## 
## Les objets suivants sont masqués depuis 'package:tseries':
## 
##     adf.test, kpss.test, pp.test
## 
## L'objet suivant est masqué depuis 'package:forecast':
## 
##     forecast
## 
## Les objets suivants sont masqués depuis 'package:fabletools':
## 
##     estimate, forecast
## 
## L'objet suivant est masqué depuis 'package:graphics':
## 
##     identify
library(forecast)
library(astsa)
library(urca)
library(seastests)
library(TSstudio)
library("cowplot")
## 
## Attachement du package : 'cowplot'
## 
## L'objet suivant est masqué depuis 'package:TSstudio':
## 
##     plot_grid
## 
## L'objet suivant est masqué depuis 'package:ggthemes':
## 
##     theme_map
## 
## L'objet suivant est masqué depuis 'package:lubridate':
## 
##     stamp
library(stats)
library(uroot)

Transformer le vecteur BD$Nombre en un objet de type ts.

La commande ts permet de convertir un vecteur en une série temporelle.

ts1 <- ts(BD$Nombre,start=c(2007,1), end=c(2016,6), frequency = 12)
ts1
##          Jan     Feb     Mar     Apr     May     Jun     Jul     Aug     Sep
## 2007  353761  325955  452936  717818 1435020 1966038 2744558 3229190 2468861
## 2008  390173  383587  464958  689467 1461643 2089183 2617998 3171273 2345464
## 2009  323236  314462  394722  680024 1309381 1892282 2512378 3003760 2266738
## 2010  274952  296528  408930  563415 1247893 1863445 2783843 3155878 2298412
## 2011  362517  333197  412673  669391 1342966 2183411 2984058 3352934 2588707
## 2012  352320  266066  360171  625837 1179714 2043692 2863095 3269202 2522887
## 2013  365604  295651  362098  557496 1467134 2372846 3263920 3885717 2963438
## 2014  404292  325952  456654  727963 1651703 2697469 4222873 4856356 3643695
## 2015  606140  509188  613091  934237 1870169 3032869 4408555 4993464 3649700
## 2016  559060  433188  628721  901585 1952667 2970996                        
##          Oct     Nov     Dec
## 2007 1473607  517954  468235
## 2008 1423374  487813  402002
## 2009 1344580  450735  410632
## 2010 1316460  458142  339591
## 2011 1428960  423107  345325
## 2012 1310986  392800  330846
## 2013 1547626  435303  402743
## 2014 1840599  668173  537728
## 2015 1852679  641458  487900
## 2016

Données représentant le nombre mensuel des visiteurs vers la Grèce de janvier 2007 à juin 2016. Chaque ligne de données représente une année soit douze valeurs Lors de la conversion du vecteur en série temporelle, on a start = c(2007,1) ,end=c(2016,6) , frequency=12 (cad mensuelle)

str(ts1)
##  Time-Series [1:114] from 2007 to 2016: 353761 325955 452936 717818 1435020 ...

Cette sortie nous indique que “ts1” est un objet de série chronologique avec 114 observations (mois).

summary(ts1)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  266066  425627  917911 1446797 2333701 4993464

La fonction de synthèse a été utilisée pour obtenir la valeur minimale, le premier quartile, la médiane, la moyenne, le troisième quartile et la valeur maximale de la série. La série va d’une valeur minimale de 1 à une valeur maximale de 4 993 464, avec une valeur médiane de 917 911. Sans contexte supplémentaire, il est difficile de fournir une interprétation plus précise de ces données.

Représenter graphiquement la série.

La représentation de la série :Ici on va diviser la fenêtre en deux parties pour bien comprendre le concept de la sérietemporelle , donc dans la première partie on va représenter le vecteur BD$Nombre , et la deuxièmepartie on va représenter la série temporelle ts1:

par (mfrow=c(1,2)) 
plot(BD$Nombre)
plot.ts(ts1)

==> La deuxième partie est la représentation de la série temporelle (les données en fonction du temps). La saisonnalité devient très apparente dans cette nouvelle série , il y a un comportement saisonnier mais pas de tendance bien évidemment puisque la variations de cette série a une amplitude qui semble augmenter au cours du temps.

Variance d’une série temporelle:

variance <- rollapply(ts1, width = 12, FUN = var, align = "right")
print(variance)
##               Jan          Feb          Mar          Apr          May
## 2007                                                                 
## 2008 1.070375e+12 1.059930e+12 1.057972e+12 1.061324e+12 1.061781e+12
## 2009 1.029116e+12 1.041304e+12 1.052582e+12 1.053655e+12 1.051369e+12
## 2010 9.410865e+11 9.441244e+11 9.419669e+11 9.549211e+11 9.543248e+11
## 2011 1.072328e+12 1.066031e+12 1.065452e+12 1.052941e+12 1.053310e+12
## 2012 1.286617e+12 1.299624e+12 1.308920e+12 1.314532e+12 1.317095e+12
## 2013 1.220537e+12 1.215079e+12 1.214751e+12 1.223479e+12 1.224538e+12
## 2014 1.716499e+12 1.709960e+12 1.691159e+12 1.664154e+12 1.665180e+12
## 2015 2.621147e+12 2.573072e+12 2.534963e+12 2.495255e+12 2.489432e+12
## 2016 2.732948e+12 2.753515e+12 2.749718e+12 2.755882e+12 2.755178e+12
##               Jun          Jul          Aug          Sep          Oct
## 2007                                                                 
## 2008 1.076729e+12 1.046322e+12 1.026861e+12 1.003022e+12 1.002007e+12
## 2009 1.026224e+12 1.001462e+12 9.459356e+11 9.308778e+11 9.289236e+11
## 2010 9.508836e+11 1.020797e+12 1.071425e+12 1.077331e+12 1.076956e+12
## 2011 1.095898e+12 1.153083e+12 1.222045e+12 1.279769e+12 1.279892e+12
## 2012 1.297319e+12 1.262143e+12 1.231766e+12 1.216850e+12 1.215387e+12
## 2013 1.277169e+12 1.401370e+12 1.645265e+12 1.749193e+12 1.747274e+12
## 2014 1.723324e+12 2.096422e+12 2.570644e+12 2.762459e+12 2.757178e+12
## 2015 2.546422e+12 2.626211e+12 2.699979e+12 2.701807e+12 2.701531e+12
## 2016 2.743448e+12                                                    
##               Nov          Dec
## 2007              1.076834e+12
## 2008 1.006562e+12 1.017338e+12
## 2009 9.341378e+11 9.328273e+11
## 2010 1.075876e+12 1.087223e+12
## 2011 1.285812e+12 1.284742e+12
## 2012 1.220279e+12 1.222794e+12
## 2013 1.738994e+12 1.724307e+12
## 2014 2.703684e+12 2.670299e+12
## 2015 2.707928e+12 2.721118e+12
## 2016
plot(variance)

Notre série montre une variance croissante au fil des années, cela peut indiquer une instabilité ou une non-stationnarité dans les données. Cela signifie que les propriétés statistiques de la série, telles que la moyenne et la variance, changent avec le temps.

ACF:

autocorr = acf(ts1, lag.max=25)

La fonction d’autocorrélation est périodique,ce qui indique une périodicité dans la série temporelle. La ligne pointillée bleue indique le niveau en-dessous duquel la corrélation n’est plus statistiquement significative. Les valeurs ≥ 0.7 sont fortement correles Les valeurs < 0.7 sont faiblement correles Les valeurs a l’interieur de l’interval de confiance ne sont pas correles. On trouve que chaque début d’année les données ont une correlation positive forte puis diminue a cause de la perte de données Les traits pointillés horizentaux représentent l’intervalle de confiance ; les valeurs à l’intérieurde cet intervalle indique que la corrélation et nulle donc les données sont indépendants , nesont pas corrélées (liées)

Vérifier si la série est corrélée:

La vérification de la corrélation entre les observations d’une série chronologique est importante car cela peut fournir des informations sur la structure de la série temporelle. Si les observations sont corrélées, cela indique qu’il existe une certaine dépendance entre les observations et que les observations passées peuvent être utilisées pour prédire les observations futures. Cela est particulièrement important pour la modélisation et la prévision de la série chronologique, car cela permet l’utilisation de modèles qui prennent en compte la dépendance temporelle des données. Pour vérifier la corrélation de notre série, on a utilisé le Box.test qui effectue un test statistique de la boîte de Ljung-Box ou si le p-value du test est inférieur à un certain seuil (généralement 0,05), on rejette l’hypothèse nulle et on conclut que la série chronologique est autocorrélée jusqu’à un certain lag. Plus précisément, le test de la boîte de Ljung-Box est utilisé pour tester l’hypothèse nulle que les corrélations croisées d’ordre k dans la série chronologique sont nulles pour tous les k jusqu’à un certain nombre de retards spécifié. Dans ce cas, le test est effectué sur les autocorrélations de la série, qui sont stockées dans la variable “autocorr$acf”.

Box.test(autocorr$acf, lag = 12, type = "Ljung-Box")
## 
##  Box-Ljung test
## 
## data:  autocorr$acf
## X-squared = 115.44, df = 12, p-value < 2.2e-16

le test de Box-Ljung a donné une statistique de test X² de 115.44 pour 12 degrés de liberté et une p-value de 2,2e-16 < 0.05 : seuil de significativité. On peut rejeter l’hypothèse nulle donc on peut conclure que les autocorrélations jusqu’à un certain lag ne sont pas nulles et que la série chronologique est autocorrélée.

Décomposition:

Commencer avec la fonction decompose(): Décomposer la série en faisant varier à chaque fois le type de la décomposition avec les options type=“add” ou “mult”, qui permettent de spécifier si on souhaite utiliser un modèle additif (par défaut) ou multiplicatif. Decompose() :Décomposer une série chronologique en composantes saisonnières, tendancielles et irrégulières en utilisant des moyennes mobiles. Traite avec un composant saisonnier additif ou multiplicatif. Le modèle additif utilisé est: Y [t] = T [t] + S [t] + e [t] Le modèle multiplicatif utilisé est: Y [t] = T [t] * S [t] * e [t] decompose() utilise la technique de la moyenne mobile pour l’estimation de tendance, d’où la présence de NA au niveau de TREND ..

Le modèle additif :

plot(decompose(ts1,"additive"))

Le modèle multiplicatif :

plot(decompose(ts1,"multiplicative"))

Analyser les résidus de chaque type de décomposition

residusadd= (decompose(ts1, type="additive"))$random 
residusmult= (decompose(ts1, type="mult"))$random
table_residus <- data.frame(residusadd, residusmult)
names(table_residus) <- c("Additive", "Multiplicative")
print(table_residus)
##        Additive Multiplicative
## 1            NA             NA
## 2            NA             NA
## 3            NA             NA
## 4            NA             NA
## 5            NA             NA
## 6            NA             NA
## 7   -299593.786      0.9448772
## 8   -319266.268      0.9559690
## 9   -173423.694      0.9696584
## 10    79115.163      1.0500113
## 11   133824.339      1.1227853
## 12   168301.102      1.2131464
## 13    72921.096      1.0724121
## 14   136202.591      1.2014925
## 15   146884.221      1.1564572
## 16   139935.716      1.1019677
## 17   157731.591      1.1107843
## 18   -40541.404      1.0149835
## 19  -402931.078      0.9171093
## 20  -344372.810      0.9621805
## 21  -255301.277      0.9503238
## 22    73041.163      1.0484153
## 23   154508.047      1.0986969
## 24   173681.602      1.0994375
## 25    90060.346      0.9469819
## 26   154847.841      1.0533496
## 27   167123.679      1.0525907
## 28   220296.924      1.1651122
## 29    96752.883      1.0681633
## 30  -148989.362      0.9848197
## 31  -423994.536      0.9401298
## 32  -430239.435      0.9715022
## 33  -258032.360      0.9748501
## 34    71188.829      1.0521937
## 35   195054.714      1.0795002
## 36   249151.352      1.1844901
## 37    85901.554      0.8343836
## 38   152009.216      1.0052782
## 39   178509.096      1.0880299
## 40    94153.924      0.9580003
## 41    21765.924      1.0070638
## 42  -189859.237      0.9605028
## 43  -166903.161      1.0297741
## 44  -300430.560      1.0026214
## 45  -250506.569      0.9695448
## 46    10082.288      1.0033400
## 47   153677.422      1.0552268
## 48   108269.227      0.9265907
## 49    92060.721      1.0319086
## 50   108368.716      1.0613634
## 51    89294.304      1.0223788
## 52    90537.091      1.0469259
## 53     3155.341      0.9938499
## 54    14992.638      1.0307736
## 55   -80927.828      1.0116906
## 56  -209215.810      0.9826466
## 57   -59384.277      1.0123886
## 58    31983.538      1.0163330
## 59    45037.630      0.9213069
## 60    70315.519      0.9114031
## 61    70711.971      0.9944370
## 62    55167.883      0.8564982
## 63    77260.096      0.9198839
## 64   111892.091      1.0276125
## 65   -85781.576      0.9231136
## 66   -49766.071      1.0207154
## 67  -127066.703      1.0268223
## 68  -223131.893      1.0097417
## 69   -61686.069      1.0349750
## 70   -23707.337      0.9773043
## 71    59268.547      0.8845462
## 72    62060.061      0.8772788
## 73    48941.971      1.0053189
## 74    -1219.242      0.8933596
## 75   -57060.779      0.8376478
## 76  -128571.159      0.8084979
## 77    11706.841      1.0012382
## 78    82823.596      1.0287994
## 79    72536.464      1.0130097
## 80   191072.940      1.0380522
## 81   172665.390      1.0489178
## 82    -7076.754      0.9860625
## 83  -123902.745      0.8355596
## 84   -87242.981      0.9155234
## 85  -156636.612      0.9423543
## 86  -253195.326      0.8198164
## 87  -269524.779      0.8713400
## 88  -277459.117      0.8677062
## 89  -133358.284      0.9224980
## 90    67251.846      0.9521886
## 91   681867.714      1.0622116
## 92   798919.565      1.0443841
## 93   481179.348      1.0339609
## 94   -89917.129      0.9395205
## 95  -269750.453      1.0279171
## 96  -332837.106      0.9814097
## 97  -303597.154      1.1604155
## 98  -351817.784      1.0970541
## 99  -332121.946      1.0398837
## 100 -250421.576      1.0123891
## 101  -71608.826      0.9615009
## 102  264451.888      0.9954290
## 103  747422.297      1.0411138
## 104  837073.649      1.0196407
## 105  404898.890      0.9921190
## 106 -144300.379      0.9135580
## 107 -347308.120      0.9611994
## 108 -411289.398      0.8774588
## 109          NA             NA
## 110          NA             NA
## 111          NA             NA
## 112          NA             NA
## 113          NA             NA
## 114          NA             NA
par (mfrow=c(1,2)) 
hist(residusadd)
hist(residusmult)

==>On remarque que les residusadd ne suivent pas la loi normale ==>Le modèle multiplicatif est meilleur quele modèle additif car les résidus du modèle multiplicatif présentent moins de fluctuations que le modèle additif. De meme l’histogramme des résidus nous conduit a dire que les données du modèle multiplicatif sont plus homogènes que celles du modèle additif.

En outre, la normalité des résidus est également importante pour l’interprétation des résultats. Si les résidus ne sont pas normalement distribués, il peut être difficile d’interpréter correctement les coefficients de régression ou les intervalles de confiance associés. En revanche, si les résidus suivent une distribution normale, cela permet une interprétation plus fiable des résultats du modèle.

Le test de Shapiro-Wilk est un test statistique qui permet de vérifier si un échantillon de données suit ou non une distribution normale. Sous R, la fonction shapiro.test() permet de réaliser ce test. Le résultat du test inclura la statistique de test W, la valeur de p et une indication sur la signification du résultat. Si la valeur de p est inférieure à 0,05, alors on peut rejeter l’hypothèse nulle selon laquelle les données suivent une distribution normale. Si la valeur de p est supérieure à 0,05, on ne peut pas rejeter l’hypothèse nulle et on peut considérer que les données suivent une distribution normale.

shapiro.test(residusmult)
## 
##  Shapiro-Wilk normality test
## 
## data:  residusmult
## W = 0.98609, p-value = 0.3656

Le résultat du test de Shapiro-Wilk indique que la statistique de test W est égale à 0,98609 et que la valeur de p est égale à 0,3656. La valeur de p étant supérieure à 0,05, on ne peut pas rejeter l’hypothèse nulle selon laquelle les résidus multiplicatifs suivent une distribution normale. Autrement dit, il n’y a pas suffisamment de preuves pour dire que les données ne suivent pas une distribution normale. Cependant, il est important de garder à l’esprit que les tests de normalité sont souvent sensibles à la taille de l’échantillon. Si l’échantillon est petit, le test peut manquer de puissance pour détecter des écarts par rapport à une distribution normale. Dans ce cas, il peut être utile d’utiliser d’autres méthodes de diagnostic pour évaluer la normalité des données.

Observation de la saisonnalité:

On va visualiser les composantes saisonnières de la série temporelle ts1 à l’aide d’un graphique en boîte (ou boxplot en anglais). La fonction “ts_seasonal” calcule les moyennes des observations de la série pour chaque période saisonnière et trace les boîtes qui montrent la variabilité des observations dans chaque période saisonnière. La visualisation permet ainsi d’observer les variations saisonnières de la série. Ensuite, On va calculer la matrice de corrélation de la série ts1. Et comme ts1 est une série temporelle univariée, la fonction renverra simplement la corrélation entre chaque observation et son observation précédente.

ts_seasonal(ts1,type="box")
ts_cor(ts1)

Graphiquement, esa boites à moustache prouve l’existence de la saisonnalité. S’il n’y avait pas d’effet saisonnier, ces 12 boites se ressembleraient ! De plus, la présence de seasonal lag 12 dans l’ACF indique la présence de la saisonnalité.

Teste de saisonnalité:

a fonction isSeasonal est une fonction R qui permet de déterminer si une série temporelle a une composante saisonnière significative. La fonction utilise une approche basée sur les filtres de Hodrick-Prescott pour décomposer la série temporelle en une tendance et une composante saisonnière. Ensuite, elle calcule une statistique de test basée sur la variance de la composante saisonnière pour déterminer si la série a une composante saisonnière significative: 1- Calcule la variance de la série temporelle originale. 2- Calcule la variance de la série temporelle différée d’une saison en utilisant la fonction diff() avec un décalage de la période de saisonnalité. 3- Compare les variances en calculant le rapport entre la variance de la série temporelle différée d’une saison et la variance de la série temporelle originale. 4- Si le rapport de variance est significativement inférieur à 1, la série est considérée comme saisonnière. Si la statistique de test est supérieure à une valeur seuil, la fonction renvoie TRUE, indiquant que la série temporelle est saisonnière. Sinon, elle renvoie FALSE, indiquant que la série n’a pas de composante saisonnière significative.

isSeasonal(ts1) #renvoie une valeur logique
## [1] TRUE

Le test de Canova et Hansen est un test de racine unitaire saisonnière. Il permet de tester si une série chronologique saisonnière suit un processus ARIMA(p,1,q)(P,1,Q)_s avec une racine unitaire. Le test est basé sur la statistique de test suivante : CH = n*(SSE1 - SSE2)/SSE2 où SSE1 et SSE2 sont les sommes des carrés des erreurs de régression pour les modèles ARIMA(p,1,q)(P,1,Q)_s avec et sans racine unitaire, respectivement, n est la taille de l’échantillon et s est la période saisonnière. La statistique de test suit une distribution chi-carré avec p+P+q+Q degrés de liberté. Sous R, le test de Canova et Hansen peut être effectué à l’aide de la fonction ch.test(). Cette fonction prend en entrée une série chronologique et retourne les résultats du test ainsi que des informations sur le modèle ARIMA estimé avec et sans racine unitaire. Voici un exemple de code :

ch.test(ts1)
## 
##  Canova and Hansen test for seasonal stability
## 
## data:  ts1
## 
##       statistic pvalue   
## Jan       0.589 0.0127 * 
## Feb      0.3277 0.1209   
## Mar      0.4101 0.0624 . 
## Apr      0.3977 0.0694 . 
## May       0.586 0.0131 * 
## Jun       0.804 0.0012 **
## Jul      0.6225 0.0089 **
## Aug      0.6079 0.0104 * 
## Sep       0.625 0.0087 **
## Oct      0.4744 0.0356 * 
## Nov      0.2825 0.1706   
## Dec      0.2333 0.2493   
## joint    1.7768 0.4102   
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
## 
## Test type: seasonal dummies 
## NW covariance matrix lag order: 12 
## First order lag: no 
## Other regressors: no  
## P-values: based on response surface regressions

Le test de Canova et Hansen pour la stabilité saisonnière évalue si les coefficients de la saisonnalité sont constants au fil du temps ou s’ils varient. La sortie montre les résultats pour chaque mois de l’année, ainsi que la statistique de test jointe pour l’ensemble des mois. Pour chaque mois, on peut voir la statistique de test et la valeur p correspondante. Si la valeur p est inférieure à un niveau de signification spécifié, cela indique que la stabilité saisonnière pour ce mois est rejetée au profit de l’hypothèse alternative de non-stationnarité saisonnière. Dans cet exemple, la stabilité saisonnière est rejetée pour les mois de janvier, mai, juin, juillet, août, septembre et octobre, car leur p-value est inférieure à 0,05. Pour les autres mois, la stabilité saisonnière n’est pas rejetée. La statistique de test jointe pour l’ensemble des mois n’est pas significative, avec une p-value de 0,4102, ce qui suggère qu’il n’y a pas de preuve de non-stationnarité saisonnière globale dans la série chronologique étudiée. Étant donné que la valeur de p est supérieure au niveau de signification de 0,05, on ne peut pas rejeter l’hypothèse nulle selon laquelle toutes les périodes ont une racine unitaire commune. Par conséquent, on peut conclure que la série chronologique a des racines unitaires saisonnières.

Plot de la composante saisonnière:

d <- decompose(ts1,"multiplicative")
d$seasonal
##            Jan       Feb       Mar       Apr       May       Jun       Jul
## 2007 0.2675508 0.2361105 0.2990129 0.4678358 0.9863914 1.5476205 2.1553155
## 2008 0.2675508 0.2361105 0.2990129 0.4678358 0.9863914 1.5476205 2.1553155
## 2009 0.2675508 0.2361105 0.2990129 0.4678358 0.9863914 1.5476205 2.1553155
## 2010 0.2675508 0.2361105 0.2990129 0.4678358 0.9863914 1.5476205 2.1553155
## 2011 0.2675508 0.2361105 0.2990129 0.4678358 0.9863914 1.5476205 2.1553155
## 2012 0.2675508 0.2361105 0.2990129 0.4678358 0.9863914 1.5476205 2.1553155
## 2013 0.2675508 0.2361105 0.2990129 0.4678358 0.9863914 1.5476205 2.1553155
## 2014 0.2675508 0.2361105 0.2990129 0.4678358 0.9863914 1.5476205 2.1553155
## 2015 0.2675508 0.2361105 0.2990129 0.4678358 0.9863914 1.5476205 2.1553155
## 2016 0.2675508 0.2361105 0.2990129 0.4678358 0.9863914 1.5476205          
##            Aug       Sep       Oct       Nov       Dec
## 2007 2.4992094 1.8797460 1.0366383 0.3407666 0.2838023
## 2008 2.4992094 1.8797460 1.0366383 0.3407666 0.2838023
## 2009 2.4992094 1.8797460 1.0366383 0.3407666 0.2838023
## 2010 2.4992094 1.8797460 1.0366383 0.3407666 0.2838023
## 2011 2.4992094 1.8797460 1.0366383 0.3407666 0.2838023
## 2012 2.4992094 1.8797460 1.0366383 0.3407666 0.2838023
## 2013 2.4992094 1.8797460 1.0366383 0.3407666 0.2838023
## 2014 2.4992094 1.8797460 1.0366383 0.3407666 0.2838023
## 2015 2.4992094 1.8797460 1.0366383 0.3407666 0.2838023
## 2016
plot(d$seasonal)

Ce plot montre que la saisonnalité varie au fil du temps ce qui nous empeche pour désaisonnaliser notre série d’utiliser la méthode de la moyenne mobile centrée qui consiste à soustraire la moyenne mobile centrée de la série chronologique et qui est simple mais elle ne fonctionne pas bien si la saisonnalité varie au fil du temps. Et de meme pour la méthode de la régression saisonnière :qui utilise une régression linéaire pour estimer les coefficients saisonniers et qui peut être utilisée que si la saisonnalité varie de manière linéaire. Par contre, on remarque que la composante saisonnière est stable car les motifs saisonniers sont cohérents d’une année à l’autre. Dans cet exemple, les valeurs saisonnières semblent être les mêmes pour chaque mois d’une année à l’autre, ce qui peut indiquer que la composante saisonnière est stable. Pour cela, nous allons effectuer la méthode STL, cette méthode utilise une décomposition en sous-séries saisonnières et tendancielles. Elle est particulièrement utile pour les séries présentant une saisonnalité non régulière ou des tendances non linéaires. Appliquons la méthode STL en utilisant la fonction stl() du package “forecast”. Cette fonction prend en entrée la série chronologique, le nombre de périodes dans une saisonnalité et le paramètre s.window qui permet de spécifier la longueur de la fenêtre de lissage pour la saisonnalité.

Désaisonnalisation:

La méthode de désaisonnalisation basée sur la décomposition STL (Seasonal and Trend decomposition using Loess) est une méthode non-paramétrique qui permet de décomposer une série temporelle en trois composantes: la tendance, la saisonnalité et la composante résiduelle. La méthode commence par déterminer une estimation de la tendance en utilisant une moyenne mobile. Ensuite, la tendance est soustraite de la série originale pour obtenir une série corrigée de tendance. La série corrigée de tendance est ensuite lissée pour extraire la saisonnalité en utilisant une méthode de lissage de Loess. Enfin, la saisonnalité est soustraite de la série corrigée de tendance pour obtenir la composante résiduelle.

stl <- stl(ts1, s.window=12)
stl
##  Call:
##  stl(x = ts1, s.window = 12)
## 
## Components
##               seasonal   trend   remainder
## Jan 2007  -988674.0806 1331039   11396.415
## Feb 2007 -1019940.0394 1333430   12464.684
## Mar 2007  -927460.3049 1335822   44574.260
## Apr 2007  -690435.5797 1338214   70039.845
## May 2007    21561.5973 1341653   71805.645
## Jun 2007   727792.8802 1345092 -106846.661
## Jul 2007  1520405.2735 1348531 -124378.078
## Aug 2007  1994107.1909 1351080 -115996.860
## Sep 2007  1167998.8493 1353629  -52766.383
## Oct 2007    67293.9154 1356177   50135.687
## Nov 2007  -894867.3186 1355721   57100.615
## Dec 2007  -972676.9043 1355264   85647.895
## Jan 2008  -993045.6335 1354807   28411.318
## Feb 2008 -1027083.0434 1349758   60912.193
## Mar 2008  -935412.6671 1344708   55662.282
## Apr 2008  -698870.7867 1339659   48678.867
## May 2008    16378.1968 1333702  111562.788
## Jun 2008   734334.6806 1327745   27103.209
## Jul 2008  1539010.1852 1321788 -242800.391
## Aug 2008  2014253.0335 1317026 -160005.848
## Sep 2008  1181023.8691 1312263 -147823.291
## Oct 2008    63170.2375 1307501   52702.732
## Nov 2008  -904542.9903 1297762   94593.588
## Dec 2008  -983459.4568 1288024   97437.683
## Jan 2009  -997574.3419 1278285   42525.197
## Feb 2009 -1034231.7325 1268946   79748.156
## Mar 2009  -943219.2440 1259606   78335.236
## Apr 2009  -706986.6099 1250266  136744.170
## May 2009    11687.7788 1245568   52124.900
## Jun 2009   741502.9300 1240870  -90091.132
## Jul 2009  1558375.0123 1236172 -282169.095
## Aug 2009  2035179.5183 1235570 -266989.638
## Sep 2009  1194850.2578 1234968 -163080.415
## Oct 2009    59755.8406 1234366   50457.965
## Nov 2009  -913601.4691 1236221  128115.069
## Dec 2009  -993809.0396 1238077  166364.434
## Jan 2010 -1007690.9168 1239932   42711.106
## Feb 2010 -1049858.1234 1243691  102695.123
## Mar 2010  -960905.8927 1247450  122385.702
## Apr 2010  -726131.8341 1251209   38337.453
## May 2010     -565.6678 1252130   -3670.915
## Jun 2010   756524.3982 1253050 -146129.182
## Jul 2010  1601383.2020 1253970  -71510.187
## Aug 2010  2080844.3353 1260403 -185369.395
## Sep 2010  1222981.7978 1266836 -191405.931
## Oct 2010    49796.2764 1273269   -6605.484
## Nov 2010  -934948.3424 1286702  106388.606
## Dec 2010 -1018246.5490 1300134   57703.285
## Jan 2011 -1018262.7516 1313567   67212.959
## Feb 2011 -1065439.4348 1326488   72148.546
## Mar 2011  -978047.1225 1339409   51311.138
## Apr 2011  -744129.8776 1352330   61190.796
## May 2011   -11070.1718 1356414   -2378.179
## Jun 2011   773837.6980 1360499   49074.683
## Jul 2011  1647226.1124 1364583  -27751.001
## Aug 2011  2129673.9174 1362577 -139317.262
## Sep 2011  1254608.1472 1360572  -26472.948
## Oct 2011    43393.4986 1358566   27000.244
## Nov 2011  -952676.4525 1348805   26978.912
## Dec 2011 -1039293.1824 1339043   45575.359
## Jan 2012 -1052559.7498 1329281   75598.643
## Feb 2012 -1110797.9999 1318868   57995.555
## Mar 2012 -1018256.3889 1308456   69971.607
## Apr 2012  -779899.1289 1298043  107693.009
## May 2012   -15097.4004 1295754 -100943.016
## Jun 2012   819748.2630 1293466  -69521.976
## Jul 2012  1715954.3855 1291177 -144036.395
## Aug 2012  2207332.7951 1297653 -235784.217
## Sep 2012  1299196.8311 1304130  -80439.665
## Oct 2012    33426.8031 1310606  -33047.050
## Nov 2012  -981203.3877 1331141   42862.850
## Dec 2012 -1072242.3274 1351675   51413.499
## Jan 2013 -1088852.3003 1372209   82247.181
## Feb 2013 -1158269.6734 1399055   54865.375
## Mar 2013 -1060696.3198 1425901   -3107.157
## Apr 2013  -818059.6905 1452748  -77191.964
## May 2013   -21676.5854 1466369   22441.543
## Jun 2013   862963.2530 1479990   29892.316
## Jul 2013  1781843.4651 1493612  -11535.284
## Aug 2013  2282060.3186 1503533  100123.357
## Sep 2013  1340762.0000 1513455  109221.170
## Oct 2013    20447.9266 1523376    3801.738
## Nov 2013 -1012731.1700 1556932 -108897.540
## Dec 2013 -1107956.5380 1590487  -79787.547
## Jan 2014 -1103717.8630 1624042 -116032.597
## Feb 2014 -1176818.0535 1671475 -168704.929
## Mar 2014 -1076576.0568 1718908 -185677.449
## Apr 2014  -831646.8927 1766340 -206730.135
## May 2014   -22162.5602 1794373 -120507.288
## Jun 2014   878688.7069 1822406   -3625.374
## Jul 2014  1806733.7213 1850438  565700.792
## Aug 2014  2310707.3362 1865702  679946.436
## Sep 2014  1358307.2550 1880966  404421.777
## Oct 2014    16493.0220 1896230  -72123.731
## Nov 2014 -1024219.0788 1911556 -219164.417
## Dec 2014 -1120927.1966 1926883 -268228.085
## Jan 2015 -1116738.5801 1942210 -219331.489
## Feb 2015 -1193880.2131 1949479 -246410.422
## Mar 2015 -1091328.1985 1956747 -252328.003
## Apr 2015  -844529.2119 1964016 -185249.556
## May 2015   -22366.3640 1966320  -73784.454
## Jun 2015   894329.5143 1968624  169915.618
## Jul 2015  1831172.5135 1970928  606454.569
## Aug 2015  2338719.8987 1959305  695438.764
## Sep 2015  1375035.0639 1947683  326982.179
## Oct 2015    11725.5715 1936060  -95106.749
## Nov 2015 -1036514.6333 1921706 -243732.880
## Dec 2015 -1134612.1812 1907351 -284838.668
## Jan 2016 -1126554.9520 1892996 -207381.233
## Feb 2016 -1206444.5114 1871952 -232319.173
## Mar 2016 -1102673.4276 1850907 -119512.756
## Apr 2016  -854581.2314 1829863  -73696.452
## May 2016   -23599.9988 1807609  168657.845
## Jun 2016   904106.7613 1785356  281533.615
seasonal_component <- stl$time.series[,1]
ts_ds <- ts1 - seasonal_component
plot(ts_ds)

isSeasonal(ts_ds)
## [1] FALSE

Test de la stationnarité:

Ce code utilise les fonctions ggAcf() et ggPacf() pour tracer la fonction d’autocorrélation (ACF) et la fonction d’autocorrélation partielle (PACF) de la série désaisonnalisée (ts_ds).

ggAcf(ts_ds) # ACF

ggPacf(ts_ds) # PACF

A première vue, on voit que celle-ci est clairement non stationnaire puisque l’ACF décroît lentement au fil du temps, ce qui indique que les observations sont fortement corrélées même à des distances temporelles et importantes et la PACF montre une décroissance lente avec des retours réguliers au-dessus du seuil de signification, ce qui indique que les observations sont fortement corrélées même à des distances temporelles importantes.

Nous pouvons en outre le confirmer avec le tests ADF : Le test Augmented Dickey-Fuller (ADF) est un test statistique utilisé pour tester la stationnarité d’une série temporelle. Il teste l’hypothèse nulle selon laquelle la série temporelle n’est pas stationnaire contre l’hypothèse alternative selon laquelle la série temporelle est stationnaire. Le test ADF peut être effectué de différentes manières, selon l’hypothèse alternative choisie. Dans ce cas, le test a été effectué avec trois hypothèses alternatives différentes : Type 1 : sans dérive ni tendance Type 2 : avec une dérive mais sans tendance Type 3 : avec une dérive et une tendance Le test ADF calcule une statistique de test appelée statistique ADF. Si cette statistique est inférieure à une certaine valeur critique, l’hypothèse nulle est rejetée, ce qui suggère que la série temporelle est stationnaire.

adf.test(ts_ds)
## Augmented Dickey-Fuller Test 
## alternative: stationary 
##  
## Type 1: no drift no trend 
##      lag     ADF p.value
## [1,]   0  0.0933   0.670
## [2,]   1 -0.3141   0.553
## [3,]   2 -0.1244   0.607
## [4,]   3  0.0805   0.666
## [5,]   4  0.1090   0.674
## Type 2: with drift no trend 
##      lag   ADF p.value
## [1,]   0 -1.94  0.3531
## [2,]   1 -3.40  0.0143
## [3,]   2 -2.61  0.0963
## [4,]   3 -1.94  0.3520
## [5,]   4 -1.88  0.3771
## Type 3: with drift and trend 
##      lag   ADF p.value
## [1,]   0 -3.07  0.1300
## [2,]   1 -5.16  0.0100
## [3,]   2 -4.32  0.0100
## [4,]   3 -3.55  0.0408
## [5,]   4 -3.59  0.0368
## ---- 
## Note: in fact, p.value = 0.01 means p.value <= 0.01

Dans les résultats affichés, les valeurs de la statistique ADF et de la p-valeur pour chaque hypothèse alternative et chaque lag sont présentées. Si la p-valeur est inférieure au niveau de signification choisi (typiquement 0,05), l’hypothèse nulle est rejetée. Dans ce cas, les résultats montrent que pour les trois types d’hypothèses alternatives, la p-valeur est supérieure au niveau de signification de 0,05. Par conséquent, on ne peut pas rejeter l’hypothèse nulle et on peut conclure que la série temporelle est non stationnaire.

Stationnariser la série:

Ce code applique la fonction diff sur la série ts_ds avec une différence d’ordre 1, c’est-à-dire qu’il calcule la différence entre chaque observation et son observation précédente. Cela permet de transformer la série en une nouvelle série, où chaque observation représente la variation entre deux observations consécutives. Cette technique est souvent utilisée pour rendre une série temporelle stationnaire, car la différenciation peut aider à éliminer la tendance et à rendre la série plus aléatoire. La nouvelle série obtenue après différenciation est stockée dans l’objet ts_ds_diff.

ts_ds_diff<- diff(ts_ds, differences=1)
autoplot(ts_ds_diff) # plot

Nous pouvons exécuter le test ADF pour vérifier la stationnarité de la nouvelle série:

adf.test(ts_ds_diff)
## Augmented Dickey-Fuller Test 
## alternative: stationary 
##  
## Type 1: no drift no trend 
##      lag   ADF p.value
## [1,]   0 -6.91    0.01
## [2,]   1 -7.45    0.01
## [3,]   2 -7.80    0.01
## [4,]   3 -6.37    0.01
## [5,]   4 -6.17    0.01
## Type 2: with drift no trend 
##      lag   ADF p.value
## [1,]   0 -6.89    0.01
## [2,]   1 -7.43    0.01
## [3,]   2 -7.78    0.01
## [4,]   3 -6.37    0.01
## [5,]   4 -6.18    0.01
## Type 3: with drift and trend 
##      lag   ADF p.value
## [1,]   0 -6.88    0.01
## [2,]   1 -7.42    0.01
## [3,]   2 -7.77    0.01
## [4,]   3 -6.36    0.01
## [5,]   4 -6.16    0.01
## ---- 
## Note: in fact, p.value = 0.01 means p.value <= 0.01

Dans ce cas, les valeurs de ADF pour tous les trois types sont inférieures à la valeur critique de 1%, ce qui signifie que nous pouvons rejeter l’hypothèse nulle selon laquelle la série n’est pas stationnaire et accepter l’hypothèse alternative selon laquelle elle est stationnaire. Les p-values étant toutes inférieures ou égales à 0,01, cela confirme également que la série est effectivement stationnaire. De plus, on peut déduire qu’il y a absence de racines unitaires.

Modélisation avec ARIMA :

ARIMA (Autoregressive Integrated Moving Average) est un modèle classique de séries temporelles qui utilise l’autorégression, la moyenne mobile et la différenciation pour modéliser une série. Il peut être utilisé pour modéliser des séries stationnaires et non-stationnaires. D’abord, nous avons divisé notre série en deux, une partie train et une partie test

train1 <- window(ts_ds_diff, end=c(2014,12))
test1 <- window(ts_ds_diff,start=c(2015,1))

La fonction auto.arima permet d’automatiquement sélectionner les ordres de l’ARIMA qui minimisent l’information de Bayes ou le critère d’information d’Akaike (AIC). Les arguments stepwise, approximation, seasonal et allowdrift sont utilisés pour spécifier les options de modélisation. stepwise=FALSE indique que la méthode stepwise de sélection de modèles ne sera pas utilisée et que l’algorithme de recherche complet sera exécuté. approximation=FALSE indique que la fonction d’approximation des erreurs ne sera pas utilisée. seasonal=FALSE indique que le modèle ne prendra pas en compte la saisonnalité. allowdrift=FALSE indique que la dérive n’est pas autorisée dans le modèle. L’argument trace=TRUE permet d’afficher les détails de l’ajustement de modèle en cours d’exécution.

ts_arima = auto.arima(train1,
                             stepwise=FALSE,
                             approximation=FALSE,
                             seasonal=FALSE,
                             allowdrift=FALSE,
                      trace=TRUE)
## 
##  ARIMA(0,0,0)            with zero mean     : 2492.408
##  ARIMA(0,0,0)            with non-zero mean : 2494.421
##  ARIMA(0,0,1)            with zero mean     : 2484.192
##  ARIMA(0,0,1)            with non-zero mean : 2486.282
##  ARIMA(0,0,2)            with zero mean     : 2485.839
##  ARIMA(0,0,2)            with non-zero mean : 2487.981
##  ARIMA(0,0,3)            with zero mean     : 2481.323
##  ARIMA(0,0,3)            with non-zero mean : 2482.657
##  ARIMA(0,0,4)            with zero mean     : 2477.914
##  ARIMA(0,0,4)            with non-zero mean : 2478.949
##  ARIMA(0,0,5)            with zero mean     : 2479.32
##  ARIMA(0,0,5)            with non-zero mean : 2480.283
##  ARIMA(1,0,0)            with zero mean     : 2485.483
##  ARIMA(1,0,0)            with non-zero mean : 2487.58
##  ARIMA(1,0,1)            with zero mean     : 2486.174
##  ARIMA(1,0,1)            with non-zero mean : 2488.314
##  ARIMA(1,0,2)            with zero mean     : 2487.151
##  ARIMA(1,0,2)            with non-zero mean : 2489.338
##  ARIMA(1,0,3)            with zero mean     : 2477.788
##  ARIMA(1,0,3)            with non-zero mean : 2478.709
##  ARIMA(1,0,4)            with zero mean     : 2479.545
##  ARIMA(1,0,4)            with non-zero mean : 2480.535
##  ARIMA(2,0,0)            with zero mean     : 2484.004
##  ARIMA(2,0,0)            with non-zero mean : 2486.119
##  ARIMA(2,0,1)            with zero mean     : 2474.829
##  ARIMA(2,0,1)            with non-zero mean : 2475.65
##  ARIMA(2,0,2)            with zero mean     : 2476.381
##  ARIMA(2,0,2)            with non-zero mean : 2477.18
##  ARIMA(2,0,3)            with zero mean     : 2478.547
##  ARIMA(2,0,3)            with non-zero mean : 2479.458
##  ARIMA(3,0,0)            with zero mean     : 2479.507
##  ARIMA(3,0,0)            with non-zero mean : 2481.501
##  ARIMA(3,0,1)            with zero mean     : 2476.289
##  ARIMA(3,0,1)            with non-zero mean : 2477.105
##  ARIMA(3,0,2)            with zero mean     : 2478.316
##  ARIMA(3,0,2)            with non-zero mean : 2479.201
##  ARIMA(4,0,0)            with zero mean     : 2480.151
##  ARIMA(4,0,0)            with non-zero mean : 2482.015
##  ARIMA(4,0,1)            with zero mean     : 2478.464
##  ARIMA(4,0,1)            with non-zero mean : 2479.357
##  ARIMA(5,0,0)            with zero mean     : 2481.017
##  ARIMA(5,0,0)            with non-zero mean : 2482.684
## 
## 
## 
##  Best model: ARIMA(2,0,1)            with zero mean
summary(ts_arima)
## Series: train1 
## ARIMA(2,0,1) with zero mean 
## 
## Coefficients:
##          ar1      ar2      ma1
##       1.0795  -0.4578  -0.8323
## s.e.  0.1072   0.0925   0.0829
## 
## sigma^2 = 1.131e+10:  log likelihood = -1233.19
## AIC=2474.38   AICc=2474.83   BIC=2484.6
## 
## Training set error measures:
##                    ME   RMSE      MAE      MPE     MAPE     MASE        ACF1
## Training set 11336.36 104645 70409.43 71.22572 205.9238 0.842597 -0.04335579

La fonction auto.arima retourne un objet Arima qui contient les coefficients de l’ARIMA, les résidus, les statistiques du modèle et autres informations. Il est possible d’accéder à ces informations en utilisant les fonctions summary, coef, fitted, resid, etc. La sortie fournit les résultats de la meilleure modélisation de la série train1 obtenue à l’aide de la fonction auto.arima. Le modèle sélectionné est un ARIMA(2,0,1) avec une moyenne nulle, ce qui indique qu’il s’agit d’un modèle autorégressif intégré à la moyenne mobile d’ordre 2 et d’ordre de différenciation 0. Les coefficients estimés pour les termes AR et MA sont également fournis, ainsi que leurs erreurs standard correspondantes. Le carré de l’écart type de l’erreur de prédiction du modèle (sigma^2) est également donné, ainsi que le logarithme de la vraisemblance du modèle et les critères d’information tels que AIC, AICc et BIC. Les mesures d’erreur sur l’ensemble d’entraînement indiquent que le modèle a une erreur moyenne (ME) de 11336.36 et un écart type racine de la moyenne (RMSE) de 104645. Le MAE est la moyenne de la valeur absolue des erreurs de prédiction et est égal à 70409.43. Le MPE est l’erreur de prédiction en pourcentage moyen et le MAPE est l’erreur de prédiction absolue en pourcentage moyen. La MASE est une mesure d’erreur de prédiction relative par rapport à un modèle de marche aléatoire, tandis que la valeur de l’ACF1 est l’autocorrélation de la première différence de la série.

Validation du modèle ARIMA(2,0,1):

Il est recommandé de toujours inspecter visuellement les résidus du modèle pour s’assurer qu’ils sont bien des bruits blancs. On peut le faire en utilisant la fonction checkresiduals. Cette fonction produit plusieurs graphiques pour vérifier si les résidus sont indépendants, ont une distribution normale et si leur variance est constante.

checkresiduals(ts_arima)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(2,0,1) with zero mean
## Q* = 10.672, df = 16, p-value = 0.8293
## 
## Model df: 3.   Total lags used: 19

Ljung-Box test est utilisé pour tester l’autocorrélation des résidus d’un modèle. Le résultat montre que la valeur de Q* est 10.672 et le degré de liberté est 16. La p-valeur est 0.8293, ce qui indique qu’il n’y a pas suffisamment de preuves pour rejeter l’hypothèse nulle selon laquelle les résidus sont indépendants et identiquement distribués. Cela suggère que le modèle ARIMA(2,0,1) avec zéro moyen est adéquat pour modéliser la série temporelle.

Effectuons maintenant les prévisions: Le modèle ARIMA(2,0,1) est utilisé pour prédire les valeurs futures de la série à l’aide de la commande predict(), en spécifiant le nombre de pas de temps à prédire avec n.ahead = length(test1). Les prévisions sont stockées dans l’objet predictions. Le premier graphe affiche les prévisions générées par le modèle, tandis que le second graphe affiche la série de données de test test1. Enfin, la fonction accuracy() est utilisée pour calculer les mesures de précision de la prédiction en comparant les prévisions avec les valeurs réelles de la série de test test1. Les mesures de précision sont renvoyées à la sortie standard.

predictions <- predict(ts_arima, n.ahead = length(test1))$pred
ggplot() + 
  geom_line(aes(x = index(test1), y = coredata(test1), color = "Test")) +
  geom_line(aes(x = index(predictions), y = predictions, color = "Forecast")) +
  scale_color_manual(values = c("Test" = "blue", "Forecast" = "darkred")) +
  labs(x = "t", y = "Yt:Nombre de touristes", color = "Legend") +
  ggtitle("Test et Forecast")

accuracy(predictions, test1)
##                ME     RMSE      MAE       MPE     MAPE      ACF1 Theil's U
## Test set 13430.73 201442.8 148301.4 -93.22257 309.4904 0.5686275 0.7097519

Le modèle ARIMA semble donner des résultats relativement précis car les valeurs de RMSE, MAE et MAPE sont assez faibles. Cependant, le MPE est assez élevé, ce qui suggère que le modèle peut avoir une tendance à sous-estimer les valeurs. Il est également important de vérifier si le coefficient d’autocorrélation est significativement différent de zéro, car cela pourrait indiquer des erreurs systématiques dans les prévisions. Le Theil’s U est inférieur à 1, ce qui indique que le modèle est meilleur que la prévision naïve.

Modélisation avec SARIMA:

SARIMA (Seasonal ARIMA) est une extension de ARIMA pour les séries saisonnières. Il utilise l’autorégression, la moyenne mobile, la différenciation saisonnière et la différenciation normale pour modéliser une série saisonnière. Commençons d’abord par stationnariser la série saisonnière ts_ds:

ts1_stat <- diff(ts1, 1)
adf.test(ts1_stat)
## Augmented Dickey-Fuller Test 
## alternative: stationary 
##  
## Type 1: no drift no trend 
##      lag   ADF p.value
## [1,]   0 -4.40    0.01
## [2,]   1 -8.08    0.01
## [3,]   2 -7.05    0.01
## [4,]   3 -5.69    0.01
## [5,]   4 -6.42    0.01
## Type 2: with drift no trend 
##      lag   ADF p.value
## [1,]   0 -4.39    0.01
## [2,]   1 -8.04    0.01
## [3,]   2 -7.02    0.01
## [4,]   3 -5.66    0.01
## [5,]   4 -6.39    0.01
## Type 3: with drift and trend 
##      lag   ADF p.value
## [1,]   0 -4.36    0.01
## [2,]   1 -8.00    0.01
## [3,]   2 -6.98    0.01
## [4,]   3 -5.62    0.01
## [5,]   4 -6.35    0.01
## ---- 
## Note: in fact, p.value = 0.01 means p.value <= 0.01

Les résultats présentés ici montrent trois types de tests ADF avec différents niveaux de tendance ou de dérive. Dans les trois types, les valeurs de la statistique ADF (colonne ADF) sont inférieures aux valeurs critiques à un niveau de signification de 0,05, ce qui suggère que la série temporelle est stationnaire. De plus, les valeurs de p (colonne p.value) sont toutes inférieures à 0,05, ce qui indique que la probabilité d’observer ces résultats par hasard est très faible.

Vérifions la saisonnalité de la nouvelle série: ce test a été défini précédemment.

isSeasonal(ts1_stat)
## [1] TRUE

Sélectionner la période d’entraînement de la série temporelle ts1_stat allant de janvier 2006 à décembre 2014 et stocker cette période dans l’objet train2. Sélectionner la période de test de la série temporelle ts1_stat allant de janvier 2015 jusqu’à la fin de la série et stocker cette période dans l’objet test2.

train2 <- window(ts1_stat, end=c(2014,12))
test2<- window(ts1_stat,start=c(2015,1))

Appliquons la fonction auto.arima() sur la série temporelle train2 pour trouver le meilleur modèle SARIMA. Les arguments stepwise=FALSE et approximation=FALSE sont utilisés pour que la fonction teste tous les modèles possibles au lieu de se limiter aux modèles les plus simples. L’argument seasonal=TRUE est utilisé pour spécifier que la série est saisonnière, et l’argument allowdrift=FALSE est utilisé pour exclure les modèles avec une dérive. L’argument trace=TRUE est utilisé pour afficher les informations de progression pendant l’exécution de la fonction. Puis, on affiche le résumé du modèle SARIMA sélectionné avec la fonction summary().

ts_sarima = auto.arima(train2,
                        stepwise=FALSE,
                        approximation=FALSE,
                        seasonal=TRUE,
                        allowdrift=FALSE,
                        trace=TRUE)
## 
##  ARIMA(0,0,0)(0,1,0)[12]                    : 2196.707
##  ARIMA(0,0,0)(0,1,1)[12]                    : 2198.693
##  ARIMA(0,0,0)(0,1,2)[12]                    : 2200.839
##  ARIMA(0,0,0)(1,1,0)[12]                    : 2198.695
##  ARIMA(0,0,0)(1,1,1)[12]                    : Inf
##  ARIMA(0,0,0)(1,1,2)[12]                    : 2203.024
##  ARIMA(0,0,0)(2,1,0)[12]                    : 2200.846
##  ARIMA(0,0,0)(2,1,1)[12]                    : Inf
##  ARIMA(0,0,0)(2,1,2)[12]                    : Inf
##  ARIMA(0,0,1)(0,1,0)[12]                    : 2197.837
##  ARIMA(0,0,1)(0,1,1)[12]                    : 2199.964
##  ARIMA(0,0,1)(0,1,2)[12]                    : 2202.104
##  ARIMA(0,0,1)(1,1,0)[12]                    : 2199.967
##  ARIMA(0,0,1)(1,1,1)[12]                    : Inf
##  ARIMA(0,0,1)(1,1,2)[12]                    : Inf
##  ARIMA(0,0,1)(2,1,0)[12]                    : 2202.05
##  ARIMA(0,0,1)(2,1,1)[12]                    : 2204.278
##  ARIMA(0,0,1)(2,1,2)[12]                    : Inf
##  ARIMA(0,0,2)(0,1,0)[12]                    : 2199.304
##  ARIMA(0,0,2)(0,1,1)[12]                    : 2201.491
##  ARIMA(0,0,2)(0,1,2)[12]                    : 2203.623
##  ARIMA(0,0,2)(1,1,0)[12]                    : 2201.494
##  ARIMA(0,0,2)(1,1,1)[12]                    : Inf
##  ARIMA(0,0,2)(1,1,2)[12]                    : Inf
##  ARIMA(0,0,2)(2,1,0)[12]                    : 2203.597
##  ARIMA(0,0,2)(2,1,1)[12]                    : Inf
##  ARIMA(0,0,3)(0,1,0)[12]                    : 2197.078
##  ARIMA(0,0,3)(0,1,1)[12]                    : 2199.304
##  ARIMA(0,0,3)(0,1,2)[12]                    : 2201.296
##  ARIMA(0,0,3)(1,1,0)[12]                    : 2199.312
##  ARIMA(0,0,3)(1,1,1)[12]                    : Inf
##  ARIMA(0,0,3)(2,1,0)[12]                    : 2201.09
##  ARIMA(0,0,4)(0,1,0)[12]                    : 2196.239
##  ARIMA(0,0,4)(0,1,1)[12]                    : 2198.316
##  ARIMA(0,0,4)(1,1,0)[12]                    : 2198.36
##  ARIMA(0,0,5)(0,1,0)[12]                    : 2196.748
##  ARIMA(1,0,0)(0,1,0)[12]                    : 2198.02
##  ARIMA(1,0,0)(0,1,1)[12]                    : 2200.167
##  ARIMA(1,0,0)(0,1,2)[12]                    : 2202.335
##  ARIMA(1,0,0)(1,1,0)[12]                    : 2200.168
##  ARIMA(1,0,0)(1,1,1)[12]                    : Inf
##  ARIMA(1,0,0)(1,1,2)[12]                    : Inf
##  ARIMA(1,0,0)(2,1,0)[12]                    : 2202.309
##  ARIMA(1,0,0)(2,1,1)[12]                    : 2204.538
##  ARIMA(1,0,0)(2,1,2)[12]                    : Inf
##  ARIMA(1,0,1)(0,1,0)[12]                    : 2199.769
##  ARIMA(1,0,1)(0,1,1)[12]                    : 2201.972
##  ARIMA(1,0,1)(0,1,2)[12]                    : 2204.146
##  ARIMA(1,0,1)(1,1,0)[12]                    : 2201.972
##  ARIMA(1,0,1)(1,1,1)[12]                    : Inf
##  ARIMA(1,0,1)(1,1,2)[12]                    : Inf
##  ARIMA(1,0,1)(2,1,0)[12]                    : 2204.09
##  ARIMA(1,0,1)(2,1,1)[12]                    : 2206.357
##  ARIMA(1,0,2)(0,1,0)[12]                    : 2193.337
##  ARIMA(1,0,2)(0,1,1)[12]                    : 2195.599
##  ARIMA(1,0,2)(0,1,2)[12]                    : 2197.705
##  ARIMA(1,0,2)(1,1,0)[12]                    : 2195.6
##  ARIMA(1,0,2)(1,1,1)[12]                    : Inf
##  ARIMA(1,0,2)(2,1,0)[12]                    : 2197.591
##  ARIMA(1,0,3)(0,1,0)[12]                    : 2195.099
##  ARIMA(1,0,3)(0,1,1)[12]                    : 2197.301
##  ARIMA(1,0,3)(1,1,0)[12]                    : 2197.322
##  ARIMA(1,0,4)(0,1,0)[12]                    : 2197.144
##  ARIMA(2,0,0)(0,1,0)[12]                    : 2199.146
##  ARIMA(2,0,0)(0,1,1)[12]                    : 2201.353
##  ARIMA(2,0,0)(0,1,2)[12]                    : 2203.46
##  ARIMA(2,0,0)(1,1,0)[12]                    : 2201.354
##  ARIMA(2,0,0)(1,1,1)[12]                    : Inf
##  ARIMA(2,0,0)(1,1,2)[12]                    : Inf
##  ARIMA(2,0,0)(2,1,0)[12]                    : 2203.375
##  ARIMA(2,0,0)(2,1,1)[12]                    : 2205.586
##  ARIMA(2,0,1)(0,1,0)[12]                    : 2192.584
##  ARIMA(2,0,1)(0,1,1)[12]                    : 2194.643
##  ARIMA(2,0,1)(0,1,2)[12]                    : 2196.703
##  ARIMA(2,0,1)(1,1,0)[12]                    : 2194.687
##  ARIMA(2,0,1)(1,1,1)[12]                    : 2196.904
##  ARIMA(2,0,1)(2,1,0)[12]                    : 2196.438
##  ARIMA(2,0,2)(0,1,0)[12]                    : 2193.279
##  ARIMA(2,0,2)(0,1,1)[12]                    : 2194.27
##  ARIMA(2,0,2)(1,1,0)[12]                    : 2194.498
##  ARIMA(2,0,3)(0,1,0)[12]                    : Inf
##  ARIMA(3,0,0)(0,1,0)[12]                    : 2200.369
##  ARIMA(3,0,0)(0,1,1)[12]                    : 2202.501
##  ARIMA(3,0,0)(0,1,2)[12]                    : 2204.704
##  ARIMA(3,0,0)(1,1,0)[12]                    : 2202.523
##  ARIMA(3,0,0)(1,1,1)[12]                    : 2204.797
##  ARIMA(3,0,0)(2,1,0)[12]                    : 2204.583
##  ARIMA(3,0,1)(0,1,0)[12]                    : 2194.479
##  ARIMA(3,0,1)(0,1,1)[12]                    : 2196.395
##  ARIMA(3,0,1)(1,1,0)[12]                    : 2196.475
##  ARIMA(3,0,2)(0,1,0)[12]                    : 2195.516
##  ARIMA(4,0,0)(0,1,0)[12]                    : 2200.986
##  ARIMA(4,0,0)(0,1,1)[12]                    : 2202.995
##  ARIMA(4,0,0)(1,1,0)[12]                    : 2203.037
##  ARIMA(4,0,1)(0,1,0)[12]                    : 2196.276
##  ARIMA(5,0,0)(0,1,0)[12]                    : 2200.956
## 
## 
## 
##  Best model: ARIMA(2,0,1)(0,1,0)[12]
summary(ts_sarima)
## Series: train2 
## ARIMA(2,0,1)(0,1,0)[12] 
## 
## Coefficients:
##          ar1      ar2      ma1
##       0.9312  -0.2398  -0.9309
## s.e.  0.1165   0.1100   0.0606
## 
## sigma^2 = 1.612e+10:  log likelihood = -1092.04
## AIC=2192.07   AICc=2192.58   BIC=2201.75
## 
## Training set error measures:
##                    ME     RMSE     MAE       MPE     MAPE      MASE        ACF1
## Training set 11056.02 116507.2 71381.5 -30.05839 78.53752 0.8132461 -0.02493748

Le meilleur modèle sélectionné est ARIMA(2,0,1)(0,1,0)[12]. Cela signifie qu’il s’agit d’un modèle autorégressif intégré de moyenne mobile saisonnier d’ordre 2, 0 et 1 pour les termes AR, I et MA, respectivement. Le modèle inclut également une différence saisonnière d’ordre 1 avec une période de 12 (saison). Les coefficients associés aux termes AR1, AR2 et MA1 sont respectivement de 0,9312, -0,2398 et -0,9309. La valeur p associée à chaque coefficient est inférieure à 0,05, ce qui suggère que les coefficients sont significatifs. La valeur de sigma² est de 1,612e+10 et la log-vraisemblance associée au modèle est de -1092,04. Les mesures d’erreur d’entraînement (ME, RMSE, MAE, MPE, MAPE, MASE et ACF1) sont également fournies, ce qui permet d’évaluer la performance du modèle sur l’ensemble d’entraînement.

Validation du modèle ARIMA(2,0,1)(0,1,0)[12]:

Validation par résidus : cette méthode consiste à vérifier si les résidus du modèle ARIMA sont aléatoires et ne présentent pas de structure temporelle. Pour cela, on peut utiliser les fonctions checkresiduals ou tsdiag, qui tracent différents graphiques pour vérifier la normalité, l’absence d’autocorrélation, l’homoscédasticité, etc.

checkresiduals(ts_sarima)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(2,0,1)(0,1,0)[12]
## Q* = 15.529, df = 16, p-value = 0.4863
## 
## Model df: 3.   Total lags used: 19

Le test de Ljung-Box est utilisé pour vérifier si les résidus d’un modèle ARIMA sont aléatoires ou s’ils présentent une autocorrélation significative. Dans ce cas, le modèle ARIMA sélectionné est (2,0,1)(0,1,0)[12]. Le résultat du test de Ljung-Box montre un Q* de 15.529 et un degré de liberté de 16, avec une p-valeur de 0.4863. Comme la p-valeur est supérieure à 0,05 (niveau de signification communément choisi), on ne peut pas rejeter l’hypothèse nulle selon laquelle les résidus sont aléatoires. Cela indique que le modèle ARIMA(2,0,1)(0,1,0)[12] est un bon ajustement pour les données.

Effectuons maintenant les prévisions:

predictions2 <- predict(ts_sarima, n.ahead = length(test2))$pred
ggplot() + 
  geom_line(aes(x = index(test2), y = coredata(test2), color = "Test")) +
  geom_line(aes(x = index(predictions2), y = predictions2, color = "Forecast")) +
  scale_color_manual(values = c("Test" = "blue", "Forecast" = "darkred")) +
  labs(x = "t", y = "Yt:Nombre de touristes", color = "Legend") +
  ggtitle("Test et Forecast")

accuracy(predictions2, test2)
##                 ME     RMSE      MAE      MPE     MAPE       ACF1 Theil's U
## Test set -11429.78 72154.85 58477.91 15.59264 24.66167 -0.2513726 0.2228339

Le tableau affiche les erreurs de prévision du modèle SARIMA pour la période de test. Les différentes métriques d’erreur sont les suivantes : Les valeurs de RMSE et MAE sont relativement élevées, ce qui indique que les erreurs de prévision ne sont pas négligeables. La valeur de MPE est positive, ce qui signifie que les prévisions sont en moyenne supérieures aux valeurs réelles. La valeur de MAPE est inférieure à 25%, ce qui indique que les prévisions ne sont pas très éloignées des valeurs réelles. La valeur de l’ACF1 est négative, ce qui suggère que les erreurs de prévision ont tendance à s’annuler mutuellement. La valeur de Theil’s U est relativement faible, ce qui indique que le modèle n’est pas très précis dans ses prévisions.

Modélisation par la méthode de Regression:

La méthode de régression est une méthode couramment utilisée pour la modélisation des séries temporelles. Elle consiste à modéliser la série temporelle en tant que fonction d’autres variables explicatives, telles que le temps, des indicateurs économiques ou météorologiques, ou d’autres facteurs pertinents. Pour appliquer la modélisation par régression linéaire sur une série temporelle, il est généralement préférable que la série temporelle soit stationnaire. La stationnarité d’une série temporelle signifie que la moyenne, la variance et la covariance ne dépendent pas du temps et restent constantes sur la période étudiée. Pour celà, nous allons appliquer la regression sur la série “ts1_stat”.

On a dévisé la série en deux sous séries train et test:

train3 <- window(ts1_stat, end=c(2014,12))
test3<- window(ts1_stat,start=c(2015,1))

La ligne de code ts_reg <- tslm(ts1_stat ~ trend + season(ts1_stat)) effectue la modélisation d’une série temporelle ts1_stat à l’aide d’un modèle de régression linéaire multiple avec une tendance linéaire et une composante saisonnière. Plus précisément, la variable dépendante est ts1_stat et les variables indépendantes sont trend et season(ts1_stat). La fonction season() est utilisée pour extraire la composante saisonnière de la série temporelle. La fonction tslm() est une fonction qui permet d’estimer les paramètres d’un modèle de régression linéaire multiple. La ligne de code suivante summary(ts_reg) affiche un résumé des résultats de la modélisation, qui incluent des informations sur les coefficients de régression, les statistiques de diagnostic, les erreurs de prévision et la qualité globale du modèle.

ts_reg <- tslm(ts1_stat ~ trend + season(ts1_stat))
summary(ts_reg)
## 
## Call:
## tslm(formula = ts1_stat ~ trend + season(ts1_stat))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -568717  -61364    8645   68317  596552 
## 
## Coefficients:
##                             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                 -28082.5    67514.0  -0.416  0.67834    
## trend                          307.5      519.0   0.592  0.55492    
## season(ts1_stat)February    -39656.5    82610.8  -0.480  0.63225    
## season(ts1_stat)March       117982.1    82596.1   1.428  0.15629    
## season(ts1_stat)April       261784.5    82584.7   3.170  0.00202 ** 
## season(ts1_stat)May         795354.9    82576.5   9.632 6.36e-16 ***
## season(ts1_stat)June        829335.8    82571.7  10.044  < 2e-16 ***
## season(ts1_stat)July        929261.6    84772.3  10.962  < 2e-16 ***
## season(ts1_stat)August      513004.5    84754.8   6.053 2.50e-08 ***
## season(ts1_stat)September  -896899.4    84740.5 -10.584  < 2e-16 ***
## season(ts1_stat)October   -1234891.3    84729.4 -14.575  < 2e-16 ***
## season(ts1_stat)November   -996793.7    84721.4 -11.766  < 2e-16 ***
## season(ts1_stat)December    -73445.3    84716.6  -0.867  0.38804    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 179700 on 100 degrees of freedom
## Multiple R-squared:  0.9436, Adjusted R-squared:  0.9368 
## F-statistic: 139.4 on 12 and 100 DF,  p-value: < 2.2e-16

Le résumé du modèle donne une liste des coefficients estimés pour chaque variable explicative incluse dans le modèle, ainsi que leurs erreurs standards, t-values et p-values. Le coefficient de détermination (R2) multiple est également fourni, qui est une mesure de la qualité de l’ajustement global du modèle. Dans ce cas, le R2 multiple est égal à 0,9436, ce qui indique que le modèle ajusté explique environ 94,36% de la variation dans les données. Le modèle indique que toutes les variables saisonnières (mois) sont significatives avec des p-values inférieures à 0,05, ce qui suggère qu’elles ont une influence significative sur la série temporelle. Cependant, le coefficient d’interception et le coefficient de tendance ne sont pas significatifs avec des p-values supérieures à 0,05. Cela signifie qu’il n’y a pas suffisamment de preuves pour conclure que la série temporelle a une tendance significative ou un niveau significatif différent de zéro.

Affichons maintenant le vecteur contenant les valeurs prédites par le modèle de régression pour chaque observation de la série temporelle ts1_stat. Il s’agit des valeurs qui minimisent l’erreur de prédiction entre les observations réelles et les prédictions du modèle.

ts_reg$fitted.values
##               Jan          Feb          Mar          Apr          May
## 2007                -67431.570    90514.530   234624.430   768502.230
## 2008   -24392.862   -63741.910    94204.190   238314.090   772191.890
## 2009   -20703.202   -60052.250    97893.850   242003.750   775881.550
## 2010   -17013.542   -56362.590   101583.510   245693.410   779571.210
## 2011   -13323.882   -52672.930   105273.170   249383.070   783260.870
## 2012    -9634.222   -48983.270   108962.830   253072.730   786950.530
## 2013    -5944.562   -45293.610   112652.490   256762.390   790640.190
## 2014    -2254.902   -41603.950   116342.150   260452.050   794329.850
## 2015     1434.758   -37914.290   120031.810   264141.710   798019.510
## 2016     5124.418   -34224.630   123721.470   267831.370   801709.170
##               Jun          Jul          Aug          Sep          Oct
## 2007   802790.630   903023.915   487074.249  -922522.196 -1260206.529
## 2008   806480.290   906713.575   490763.909  -918832.536 -1256516.869
## 2009   810169.950   910403.235   494453.569  -915142.876 -1252827.209
## 2010   813859.610   914092.895   498143.229  -911453.216 -1249137.549
## 2011   817549.270   917782.556   501832.889  -907763.556 -1245447.889
## 2012   821238.930   921472.216   505522.549  -904073.895 -1241758.229
## 2013   824928.590   925161.876   509212.209  -900384.235 -1238068.569
## 2014   828618.250   928851.536   512901.869  -896694.575 -1234378.909
## 2015   832307.910   932541.196   516591.529  -893004.915 -1230689.249
## 2016   835997.570                                                    
##               Nov          Dec
## 2007 -1021801.529   -98145.640
## 2008 -1018111.869   -94455.980
## 2009 -1014422.209   -90766.320
## 2010 -1010732.549   -87076.660
## 2011 -1007042.889   -83387.000
## 2012 -1003353.229   -79697.340
## 2013  -999663.569   -76007.680
## 2014  -995973.909   -72318.020
## 2015  -992284.249   -68628.360
## 2016
autoplot(ts_reg$fitted.values)

Validation de la méthode de regression:

checkresiduals(ts_reg)

## 
##  Breusch-Godfrey test for serial correlation of order up to 23
## 
## data:  Residuals from Linear regression model
## LM test = 78.478, df = 23, p-value = 5.59e-08

Le test de Breusch-Godfrey est utilisé pour détecter la présence de corrélation série temporelle dans les résidus d’un modèle de régression. Dans ce cas, le test a été effectué sur les résidus d’un modèle de régression linéaire. La statistique de test LM est de 78,478 avec 23 degrés de liberté, et la valeur de p est de 5,59e-08. Puisque la valeur de p est inférieure à 0,05, on rejette l’hypothèse nulle selon laquelle il n’y a pas de corrélation série temporelle dans les résidus, et on conclut qu’il y a une corrélation série temporelle dans les résidus. Cela suggère que le modèle de régression linéaire n’est pas suffisant pour modéliser la série temporelle.

Modélisation par la méthode de lissage exponentielle:

Les modèles de lissage exponentiel sont des méthodes de prévision de séries temporelles qui s’appuient sur la notion de moyenne mobile exponentielle pour estimer les tendances locales et saisonnières d’une série temporelle. Ces modèles ne sont pas linéaires, car ils utilisent des coefficients de lissage exponentiels pour donner plus ou moins d’importance aux observations passées en fonction de leur éloignement dans le temps. Les modèles de lissage exponentiel peuvent être appliqués dans des situations où la série temporelle a une tendance, une saisonnalité et des composantes d’erreur qui sont stables ou qui changent lentement au fil du temps. En d’autres termes, ils sont adaptés aux séries temporelles qui ne présentent pas de tendances ou de motifs de saisonnalité complexes et qui ne présentent pas de changements brusques dans les modèles de comportement. Le lissage exponentiel est souvent utilisé pour modéliser des séries temporelles non stationnaires et saisonnières. La série doit être lissée pour éliminer les variations de courte durée et mettre en évidence les tendances sous-jacentes qui peuvent être utilisées pour la prévision. Les modèles de lissage exponentiel sont donc appliqués directement sur la série temporelle originale.

On divise d’abord la série originale en deux sous séries:

train4 <- window(ts1, end=c(2014,12))
test4 <- window(ts1,start=c(2015,1))

La commande ts1_hw <- HoltWinters(train4, beta=TRUE, gamma=TRUE) applique le modèle de Holt-Winters avec tendance et saisonnalité additive sur la série train4 pour estimer les paramètres nécessaires à la prévision des valeurs futures. La tendance est modélisée par l’exponentielle double, et la saisonnalité est prise en compte par la méthode de lissage exponentiel triple. Les paramètres beta et gamma permettent de déterminer si la tendance et/ou la saisonnalité sont modélisées de manière additive ou multiplicative.

ts1_hw <- HoltWinters(train4, beta=TRUE, gamma=TRUE)
predictions4 <- predict(ts1_hw, n.ahead=18)
ggplot() + 
  geom_line(aes(x = index(test4), y = coredata(test4), color = "Test")) +
  geom_line(aes(x = index(predictions4), y = predictions4, color = "Forecast")) +
  scale_color_manual(values = c("Test" = "blue", "Forecast" = "darkred")) +
  labs(x = "t", y = "Yt:Nombre de touristes", color = "Legend") +
  ggtitle("Test et Forecast")

accuracy(predictions4, test4)
##               ME    RMSE     MAE      MPE     MAPE      ACF1 Theil's U
## Test set 1005335 1170533 1005335 78.99331 78.99331 0.7117418  2.075206

ME : Erreur moyenne, c’est la moyenne des différences entre les valeurs observées et prédites. Ici, elle vaut 1005335, ce qui signifie que les prédictions ont tendance à sous-estimer les valeurs observées. RMSE : Erreur quadratique moyenne, c’est la racine carrée de la moyenne des carrés des différences entre les valeurs observées et prédites. Ici, elle vaut 1170533, ce qui représente la déviation standard des erreurs entre les valeurs observées et prédites. MAE : Erreur absolue moyenne, c’est la moyenne des valeurs absolues des différences entre les valeurs observées et prédites. Ici, elle vaut 1005335, ce qui indique une tendance de sous-estimation des prédictions. MPE : Erreur de pourcentage moyen, c’est la moyenne des différences en pourcentage entre les valeurs observées et prédites. Ici, elle vaut 78.99331, ce qui indique que les prédictions sont en moyenne inférieures de 78,99% aux valeurs observées. MAPE : Erreur absolue de pourcentage moyen, c’est la moyenne des valeurs absolues des différences en pourcentage entre les valeurs observées et prédites. Ici, elle vaut également 78.99331, ce qui confirme la tendance de sous-estimation des prédictions. ACF1 : Le coefficient de corrélation autocorrélatif de la première différence. S’il est proche de 1 ou -1, cela indique une forte autocorrélation dans les données. Ici, il vaut 0.7117418, ce qui indique une certaine autocorrélation dans les résidus du modèle. Theil’s U : Coefficient d’efficacité de la prédiction, il mesure l’efficacité de la prédiction par rapport à une méthode de référence telle que la prédiction naïve. Ici, il vaut 2.075206, ce qui indique que le modèle de lissage exponentiel est deux fois plus efficace que la méthode de prédiction naïve.

Comparaison des 3 modèles:

ARIMA_res = data.frame(accuracy(predictions, test1))
SARIMA_res = data.frame(accuracy(predictions2, test2))
LISSAGE_res = data.frame(accuracy(predictions4, test4))
resf <- data.frame(MAE=c(ARIMA_res$MAE,SARIMA_res$MAE,LISSAGE_res$MAE),MAPE = c(ARIMA_res$MAPE, SARIMA_res$MAPE, LISSAGE_res$MAPE),RMSE=c(ARIMA_res$RMSE,SARIMA_res$RMSE,LISSAGE_res$RMSE))
rownames(resf) <- c("ARIMA", "SARIMA", "LISSAGE EXPONENTIELLE")
resf
##                              MAE      MAPE       RMSE
## ARIMA                  148301.44 309.49035  201442.75
## SARIMA                  58477.91  24.66167   72154.85
## LISSAGE EXPONENTIELLE 1005335.43  78.99331 1170533.25

==> Dans ce cas, le modèle SARIMA a les valeurs les plus faibles pour les trois mesures, ce qui indique qu’il a la meilleure performance de prévision parmi les trois modèles testés. Le modèle de lissage exponentiel a les valeurs les plus élevées pour les trois mesures, ce qui indique qu’il a la pire performance de prévision.

PDF <- predict(ts_sarima, n.ahead=48)
orig_series <- cumsum(PDF$pred)
orig_series1 <- ts(orig_series,start=c(2015,1), end=c(2018,12), frequency = 12)
autoplot(orig_series1)